home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / graph_alg / _mcmflow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  8.9 KB  |  337 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  _mcmflow.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16.  
  17. //------------------------------------------------------------------------------
  18. //
  19. // MIN_COST_MAX_FLOW(G,s,t,cap,cost,flow)
  20. //
  21. // computes a maximal s-t flow of minimal total cost in network (G,cap,cost)
  22. //
  23. // Algorithm: successive shortest path augmentation + capacity scaling
  24. //
  25. // Based on:  J. Edmonds and R. M. Karp: 
  26. //            "Theoretical Improvements in Algorithmic Efficiency for Network 
  27. //            Flow Problems", Journal of the  ACM, Vol. 19, No. 2, 1972
  28. //
  29. //            R.K. Ahuja, T.L. Magnanti, J.B. Orlin: 
  30. //            "Network Flows", Section 10.2, Prentice Hall Publ. Company, 1993
  31. //
  32. //
  33. // S. N"aher  (December 1993)
  34. //------------------------------------------------------------------------------
  35.  
  36. #include <LEDA/graph_alg.h>
  37. #include <LEDA/node_pq.h>
  38.  
  39.  
  40. static node DIJKSTRA_IN_RESIDUAL(const graph& G, node s, int delta, 
  41.                                  const edge_array<num_type>& cost,
  42.                                  const edge_array<num_type>& cap_minus_flow,
  43.                                  const edge_array<num_type>& flow,
  44.                                  const node_array<num_type>& excess,
  45.                                  node_array<num_type>& pi,
  46.                                  node_array<edge>& pred)
  47. {
  48.    /*  Computes a shortest path from node s until a node t with excess <= -delta
  49.        is reached in the delta-residual graph of network G. Edges are traversed
  50.        in both directions using the forall_in/out_edges iteration macros.The 
  51.        residual capacity of an edge e is cap_minus_flow[e] if e is used from 
  52.        source(e) to target(e) and flow[e] otherwise. The reduced cost of an
  53.        edge e = (u,v) is equal to cost[e] + pi[u] - pi[v] if it is used from 
  54.        u and -cost[e] + pi[v] - pi[u] if it is used from v.
  55.  
  56.        precondition: every edge in the delta-residual network has non-negative 
  57.                      reduced cost.
  58.        
  59.        returns t (nil,if there is no path from s to a node with excess<-delta)
  60.  
  61.        and
  62.  
  63.        pred[v] = predecessor edge of v on the shortest path from s to t
  64.  
  65.        pi[v]   = pi'[v] + dist[v] - dist[t], if v is permanently labeled
  66.                  pi'[v],                     if v is temporarily labeled
  67.                  (here pi' is the old potential vector)
  68.  
  69.   */
  70.  
  71.   node_array<num_type> dist(G,MAXINT);
  72.   node_pq<num_type> PQ(G);
  73.   node v;
  74.   edge e;
  75.  
  76.   node t = nil;
  77.  
  78.   node_list perm_labeled; // list of permanently labeled nodes
  79.  
  80.   dist[s] = 0;
  81.   PQ.insert(s,0);
  82.  
  83.   while (!PQ.empty())
  84.   { node u = PQ.del_min();
  85.  
  86.     perm_labeled.append(u);
  87.  
  88.     if (excess[u] <= -delta) 
  89.     { t = u;
  90.       break;
  91.      }
  92.  
  93.     num_type du = dist[u] + pi[u]; 
  94.  
  95.     forall_out_edges(e,u) 
  96.       if (cap_minus_flow[e] >= delta)  // e in delta-residual graph
  97.       { v = target(e);
  98.         num_type c = du - pi[v] + cost[e];
  99.         if (c < dist[v])
  100.         { if (dist[v] == MAXINT) // v has not been reached before
  101.             PQ.insert(v,c);
  102.           else
  103.             PQ.decrease_inf(v,c);
  104.           dist[v] = c;
  105.           pred[v] = e;
  106.          }
  107.        }
  108.  
  109.     forall_in_edges(e,u) 
  110.       if (flow[e] >= delta)  // e in delta-residual graph
  111.       { v = source(e);
  112.         num_type c = du - pi[v] - cost[e];
  113.         if (c < dist[v])
  114.         { if (dist[v] == MAXINT) // v has not been reached before
  115.             PQ.insert(v,c);
  116.           else
  117.             PQ.decrease_inf(v,c);
  118.           dist[v] = c;
  119.           pred[v] = e;
  120.          }
  121.        }
  122.  
  123.    }
  124.  
  125.    if (t != nil) // update potentials of all permanently labeled nodes
  126.    { num_type dt = dist[t];
  127.      forall(v,perm_labeled) pi[v] += dist[v] - dt;
  128.     }
  129.  
  130.    return t;
  131. }
  132.  
  133.  
  134.  
  135. num_type MIN_COST_MAX_FLOW(graph& G, node src, node sink,
  136.                                               const edge_array<num_type>& cap0,
  137.                                               const edge_array<num_type>& cost0,
  138.                                               edge_array<num_type>& flow0)
  139. { 
  140.   double n = G.number_of_nodes();
  141.   double m = G.number_of_edges();
  142.  
  143.   num_type total_flow = 0;               // total flow value (returned)
  144.   int  count = 0;                        // number of augmentations
  145.   node v;
  146.   edge e;
  147.  
  148.   num_type supply = MAX_FLOW(G,src,sink,cap0,flow0);
  149.  
  150.  
  151.   num_type C = 0;  // maximal cost of an s-t path in G
  152.  
  153.   forall_edges(e,G)
  154.   { if (cost0[e] > C)  C = cost0[e];
  155.     if (cost0[e] < -C) C = -cost0[e];
  156.    }
  157.  
  158.   C = C * G.number_of_nodes();
  159.  
  160.  
  161.  
  162.   list<edge> art_edges;  // list of artificial edges (see below)
  163.   list<edge> neg_edges;  // list of negative cost edges
  164.  
  165.   // add artificial edges to guarantee existence of a path of infinite
  166.   // capacity (MAXINT) and huge cost (C) between any pair of nodes  (v,w)
  167.   // with b(v) > 0 and b(w) < 0
  168.    
  169.   //forall_nodes(v,G) 
  170.   //if (v != src) 
  171.   //{ art_edges.append(G.new_edge(v,src));
  172.   //  art_edges.append(G.new_edge(src,v));
  173.   // }
  174.  
  175.   art_edges.append(G.new_edge(src,sink));
  176.  
  177.  
  178.   node_array<num_type> pi(G,0);           // node potential
  179.   node_array<num_type> excess(G,0);       // excess
  180.   node_array<edge>     pred(G,nil);       // predecessor edge on shortest path
  181.   edge_array<num_type> flow(G,0);
  182.   edge_array<num_type> cap(G,0);
  183.   edge_array<num_type> cost(G,0);
  184.   edge_array<num_type> cap_minus_flow(G,0);
  185.  
  186.  
  187.   forall_edges(e,G)
  188.   { cap[e] = cap0[e];
  189.     cost[e] = cost0[e];
  190.     cap_minus_flow[e] = cap[e];
  191.     if (cost[e] < 0) neg_edges.append(e);
  192.    }
  193.  
  194.   forall(e,art_edges)
  195.   { cap[e] = MAXINT;
  196.     cost[e] = C;
  197.     cap_minus_flow[e] = MAXINT;
  198.    }
  199.  
  200.   excess[src]  =  supply;
  201.   excess[sink] = -supply;
  202.  
  203.  
  204.   // eliminate negative cost edges by edge reversals 
  205.   // (Ahuja/Magnanti/Orlin, section 2.4)
  206.  
  207.   forall(e,neg_edges)
  208.   { node u = source(e);
  209.     node v = target(e);
  210.     excess[u] -= cap[e];
  211.     excess[v] += cap[e];
  212.     cost[e] = -cost[e];
  213.     G.rev_edge(e);
  214.   }
  215.  
  216.   int delta = 1;
  217.  
  218.   double delta_max = supply * m/(n*n);  // seems to be a good choice
  219.  
  220.   while (2*delta <= delta_max) delta = 2*delta;
  221.  
  222.   while (delta > 0)  // delta scaling phase
  223.   {
  224.     forall_edges(e,G)  
  225.     { // Saturate all edges entering the delta residual network which have
  226.       // a negative reduced edge cost. Then only the reverse edge (with positive
  227.       // reduced edge cost) will stay in the residual network. 
  228.  
  229.       node u = source(e);
  230.       node v = target(e);
  231.       num_type c = cost[e] + pi[u] - pi[v];
  232.  
  233.       if (cap_minus_flow[e] >= delta && c < 0)
  234.         { excess[v] += cap_minus_flow[e];
  235.           excess[u] -= cap_minus_flow[e];
  236.           flow[e] = cap[e];
  237.           cap_minus_flow[e] = 0;
  238.          }
  239.       else
  240.         if (flow[e] >= delta && c > 0)
  241.         { excess[u] += flow[e];
  242.           excess[v] -= flow[e];
  243.           flow[e] = 0;
  244.           cap_minus_flow[e] = cap[e];
  245.          }
  246.  
  247.      }
  248.  
  249.  
  250.     // As long as there is a node s with excess >= delta compute a minimum
  251.     // cost augmenting path from s to a sink node t with excess <= -delta in 
  252.     // the delta-residual network and augment the flow by at least delta units 
  253.     // along P (if there are nodes with excess > delta and excess < -delta
  254.     // respectively, the existence of P is guaranteed by the artificial edges 
  255.     // inserted at the begining of the algorithm)
  256.  
  257.     node s;
  258.     node t;
  259.   
  260.     forall_nodes(s,G)
  261.     while (excess[s] >= delta)
  262.     { 
  263.       count++;
  264.   
  265.       t = DIJKSTRA_IN_RESIDUAL(G,s,delta,cost,cap_minus_flow,flow,excess,pi,pred);
  266.       if (t == nil) break; // there is no node with excess < -delta
  267.   
  268.       // compute maximal amount f of flow that can be pushed through P
  269.   
  270.       num_type f = excess[s];
  271.   
  272.       v = t; 
  273.       while (v != s)
  274.       { e = pred[v]; 
  275.         num_type rcap;
  276.         if (v==target(e))
  277.            { rcap = cap_minus_flow[e];
  278.              v = source(e);
  279.             }
  280.         else
  281.            { rcap = flow[e];
  282.              v = target(e);
  283.             }
  284.         if (rcap < f) f = rcap;
  285.        }
  286.   
  287.       if (f > -excess[t]) f = -excess[t];
  288.   
  289.   
  290.       // add f units of flow along the augmenting path 
  291.   
  292.       v = t; 
  293.       while (v != s)
  294.       { e = pred[v]; 
  295.         if (v==target(e))
  296.            { flow[e] += f;
  297.              cap_minus_flow[e] -= f;
  298.              v = source(e);
  299.             }
  300.         else
  301.            { flow[e] -= f;
  302.              cap_minus_flow[e] += f;
  303.              v = target(e);
  304.             }
  305.        }
  306.   
  307.        excess[s] -= f;
  308.        excess[t] += f;
  309.   
  310.     } // end of delta-phase
  311.   
  312.    delta /= 2;
  313.   
  314.   } // end of scaling
  315.   
  316.  
  317.  
  318.   // restore negative edges
  319.   forall(e,neg_edges) 
  320.   { flow[e] = cap[e] - flow[e];
  321.     G.rev_edge(e);
  322.    }
  323.  
  324.   // delete artificial edges
  325.   forall(e,art_edges) G.del_edge(e); 
  326.  
  327.   // total flow  = total flow out of the source node
  328.   forall_out_edges(e,src) total_flow += flow[e];
  329.  
  330.   forall_edges(e,G) flow0[e] = flow[e];
  331.  
  332.   //cout << "count = " << count << endl;
  333.  
  334.   return total_flow;
  335. }
  336.  
  337.